% This script generates ensembles of temperature trajectories for figure 2. The same ensembles are used for table 2 in the SI.
% This script also generates ensembles of temperature trajectories for figure 3 in the SI. See options below.
close all
clear all
 
% Load RCP data (This script is run once for each RCP. Each time the script is run, you must "comment out" all but one of the RCPs below.)
%load rf_rcp8_5.txt
%load rf_rcp6_0.txt
%load rf_rcp4_5.txt
load rf_rcp2_6.txt
%yrs=rf_rcp8_5(:,1);
%yrs=rf_rcp6_0(:,1);
%yrs=rf_rcp4_5(:,1);
yrs=rf_rcp2_6(:,1);
%forcing=rf_rcp8_5(:,2);
%forcing=rf_rcp6_0(:,2);
%forcing=rf_rcp4_5(:,2);
forcing=rf_rcp2_6(:,2);

% Set start time in 2020
tstart=256;

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Specify initial temperature anomaly in 2020
DeltaT0=1; % degrees C.

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Specify model parameters %%%%%%%%%%%%%%%%%%
% We will measure time in years,
% Specify size of time interval between points.
% so need a factor of number of seconds in year
% Effectively divides time in seconds by yearlength
dt=1; 
% So then have to multiply speed by yearlength
yearlength=3.1 *1e7;

% Specify number of points at which solution is wanted
nPeriods=length(forcing(tstart:end)); % Same as RCP dataset

% Specify # of intermediate steps taken between the sample points
nSteps=12; % Monthly

% Number of temperature trajectories without noise
Ntrials0=1;

% Number of temperature trajectories with noise
Ntrials=10000;

% Effective heat capacity (Figures 2 and 3 in the SI plot trajectories that assume different heat capacities. Comment out lines below as needed.)
%C=0.8*1e9; % Watts/metre^2 /Kelvin. 
C=1*1e9; % Watts/metre^2 /Kelvin 
%C=1.2*1e9; % Watts/metre^2 /Kelvin. 

% Standard deviation of noise (Figures 2 and 3 in the SI plot trajectories that assume different amounts of noise. Comment out lines below as needed.)
% sigmaQ = 0.6725*1e8;
sigmaQ = 0.9375*1e8;
% sigmaQ = 1.0625*1e8;

% Equilibrium climate sensitivity and feedback parameter
CS = 3; % degrees C
lambda=3.7 / CS; % Joules/metres^2/Kelvin

% Functional parameters
Speed = lambda*yearlength/C; % Kelvin per year 
Level=ts2func(forcing(tstart:end)/lambda); % Convert F/lambda to a function

% Generate trajectories without noise
Sigma=0;
[DeltaTnoiseless,Times0]=Hasselmann(Speed,Level,Sigma,DeltaT0,Ntrials0,nPeriods,nSteps,dt);

% Generate trajectories with noise
Sigma=sigmaQ/C;
[DeltaT,Times]=Hasselmann(Speed,Level,Sigma,DeltaT0,Ntrials,nPeriods,nSteps,dt);

%Create a CSV write file
% First suppress the unit 3rd dimension in Delta T matrix
Temp=squeeze(DeltaT);
% Then concatenate 
M=[Times DeltaTnoiseless Temp];
 
% filename = ['RCP26output_lowSigma.txt']; % Write output to this file when assuming C=0.8*1e9 and sigmaQ = 0.6725*1e8. Change "26" to "45", "60", or "85" when running other RCPs.
filename = ['RCP26output.txt'];  % Write output to this file when assuming C=1*1e9 and sigmaQ = 0.9375*1e8. Change "26" to "45", "60", or "85" when running other RCPs.
% filename = ['RCP26output_highSigma.txt']; % Write output to this file when assuming C=1.2*1e9 and sigmaQ = 1.0625*1e8. Change "26" to "45", "60", or "85" when running other RCPs.
csvwrite(filename,M)